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Abstract 



We briefly describe the Generalized Schur Algorithm for the superfast solution of 
positive definite Toeplitz systems of equations and its relationship with Schur’s 
algorithm and the Szego recursions. We then present some expeiimental results 
obtained with our FORTRAN implementation of this superfast Toeplitz solver. 
We will see that the algorithm displays favorable behavior in that the growth 
rates of the resulting residuals and errors are comparable with those of the Szego 
recursions. 



1. Introduction 

The problem of efliciently solving a linear system of equations M x = b where 
M is an Hermitian positive definite Toeplitz matrix is of considerable interest in 
many areas of mathematics and engineering. While standard Gauss or Choleski 
methods for solving such a system of order n require O(ii^) arithmetic operations, 
there are several algorithms for obtaining the solution in 0(n‘) operations 
[21,24,5,18,13]. We refer to these algorithms as fast Toeplitz solvers. More 
recently, algorithms have been presented for solving a Toeplitz system of equa- 
tions using 0{n\og^n) arithmetic operations [11,6,7,15,17,22]. These siiperfast 

t Presented at the International Conference on Linear Algebra and Its Ai)plications, 
Valencia, Spain, Sej^teinber, 1987. d'his research was su|)portcd in part by the National 
Science Foundation under grant DMS-8704196. 
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Toeplilz solvers are all based on certain doubling procedures. While they are moi'e 
complicated to implement than the fast methods, they require substantially less 
computational work for sufficiently large problems. 

Current fast and superfast Toeplitz solvers are numerically unstable, and 
therefore unreliable, when applied to arbitrary Toeplitz systems [8]. However, 
Cybenko [9, 10] has shown that the algorithms of Levinson, Durbin, and Trench 
are numerically stable for the class of positive definite Toeplitz matrices in the 
sense that accurate results can be expected when the problem is well conditioned. 
Stability of some superfast algorithms can also be expected for the class of posi- 
tive definite matrices [8]. We remark that some fast and superfast algorithms are 
also potentially forward stable for other classes of matrices, for example, for 
totally positive Toeplitz matrices. 

The numerical reliability of fast and superfast Toeplitz solvers is a topic of 
ongoing research. Since they are relatively neAv on the scene and, moreover, 
implementations of these algorithms are less tractable than those of the fast algo- 
rithms, numerical results for superfast Toeplitz solvers are currently rare. 
Ivnowledge of the numerical reliability is important because superfast methods 
allow for the solution of problems which would require a prohibitive amount of 
time for a fast method. We are not aware of any experimental results for any 
superfast Toeplitz solvers except for those presented in [2]. 

Positive definite Toeplitz matrices play an important role in several aspects 
of classical analysis including the classical theory of Szego polynomials, the tri- 
gonometric moment problem, and the Carathebdory coefficient problem [l, 16]. 
Several Toeplitz solvers are direct manifestations of algorithms and formulas that 
were developed in the classical theory (see, e.g., [20, 19, 3]). In [3], the superfast 
Toeplitz solver presented in [17] and [22] is described in terms of a generalization 
of Schur’s classical algorithm. 

In this paper we give a general description of the Generalized Schur Algo- 
rithm and its use in the superfast solution of Toeplitz systems. We then present 
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some results of our experimental code for the superfast solution of real positive 
definite Toeplitz matrices. Refer to [3] for a development of the Generalized Schur 
Algorithm and its relation with fast Toeplitz solvers, and see [2,-l] for descriptions 
of the implementation of the algorithm for complex and real positive definite Toe- 
plitz matrices. 



2. Some Toeplitz Solvers 

The Toeplitz solvers we consider can be divided into two phases: a factorization 
phase and a solution phase. In this section we summarize some fast Toeplitz 
solvers based on the Szego recursions and Schur’s algorithm. We then present the 
generalization of Schur’s algorithm for the superfast solution of a positive definite 
Toeplitz system of equations. 

2.1. The Szego Recursions. The Hermitian positive definite Toeplitz matrix 
M= ~ 1/^y-t- ]/t= 0 ~ defines an inner product on the set of polynomi- 

als of degree at most by 

The monic polynomials Xy, degXy =j, j = Q, ' ' d* fli^t are orthogonal with 

respect to this inner product are the (monic) Szego polynomials determined by A/. 
These polynomials can be constructed using roughly multiplications and addi- 
tions by the Szego recurrence relations. 
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Szego Recursions: 



%0 • ~ ~ /"''O’ 

for j = 0, 1, . . . , ?i —1 

7;+i : = -<l,Xxy>/<5y 

Xy+i(X) : = Xxy(X) + 7y+i W 

^y+i: = ^y(l-hy^iP) 



where Xy(X) : = X-^ x(l/X) is the polynomial obtained by conjugating and revers- 
ing the coefficients of Xy- The quantities 7y are called Schur parameters', they 
serve as intermediate quantities in the Szego recursions and are often of physical 
and mathematical significance. 

The Szego recursions comprise the first phase of several Toeplitz solvers. In 
fact, they are equivalent with the Levinson-Durbin for the solution of the X^ule- 
Walker equations of linear prediction. Define R to be the unit right triangular 
matrix whose kih column contains the coefficients of Xk (0^^ ^^0- 
definition of the Szego polynomials we have 

R*MR=D, M-'^ = R D-'^ R-* , 

where D := diag[6^]o” and df. := <Xjt)X7>- Thus, the Szego polynomials deter- 
mine the reverse Choleski factorization of This is the factorization that 

underlies the solution phase of Levinson’s algorithm for the solution of a Toeplitz 
system of equations. 

In Trench’s algorithm a relationship between the entries of M~^ is used to 
construct M~^ from its last column, which is given by Xn ^n- This relation- 
ship follows from the Christoffel-Darboux formula for Szego polynomials [24, 20]. 
The matrix interpretation of the Christoffel-Darboux-Szego formula is the 
Gohberg-Semencul formula 
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K = 'i\ r{ - To To 

where 

0 0 ■ 0 

Po ^ 

; • . • . 0 

Pn-l ■ ■ ■ Po 0 

1 Pn-l • • • Po ’ 

0 1 '•.••• 

Pn — 1 

0 • • • 0 1 

and Xnl^) = YjO Pj ’ (Pn = l)- See, e.g., [12]. Thus, the Gohberg-Semencul for- 
mula provides a factorization of that is determined by the last Szego polyno- 
mial Xn squared norm <5„ . The product M~^b can then be obtained in 

O(7ilog?i) operations by using Fast Fourier Transform (FFT) techniques. This 
provides an 0{n logji ) procedure for the second phase of a Toeplitz solver. 



To 

T, 



2.2. Schur’s Algorithm. Another cla.ss of Toeplitz solvere is based on the clas- 
sical algorithm of Schur [23]. Let be a Schur function, which is a holomorphic 
mapping of the open unit disk in the complex plane into its closure. Schur’s algo- 
rithm generates a continued fraction representation of a given Schur function 
(f> = 4>o by the successive application of linear fractional transformations (LFTs). 
Given a Schur function </>„_!, define 7,, = and 






^n-li>')-ln 

1 - 



It is easy to see that (/>„ is also a Schur function. Furthermore, if [7,; |=1 then 
in this case </> is a rational function and the algorithm terminates. 
Otherwise, we repeat the procedure to obtain + i from . In this way we 
obtain the Schur continued fraction representation of S, 
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<f> = Tn[4>n) = 



tn{K)- 



where 



kij) = 



In 

1 + X r 



The function T„ (0) is referred to as the nth approximant of (f), and is called the 
7ith tail of 4>. Thus, Schur’s algorithm effects a sequence of elementary linear frac- 
tional transformations on the initial Schur function <p, and the Schur parameters 
7y are intermediate quantities that determine the elementary linear fractional 
transformations. 

To implement Schur’s algorithm we write each Schur function as the quotient 

^n X) ' 

of formal power series, (/»„ = = — with (3^ q> 0 as a, partial nor- 

Pn Yjj^o^nj' 

malization. Organizing Schur’s algorithm so that the coefficient pairs Qq k 
entered and processed sequentially leads to the progressive form of Schur algo- 
rithm, which is given below. 



Schur’s Algorithm: 

for ^ = 1, 2, 3, • •• while \^f. \ < 1 
enter /?o,a-i 

for y = 1, 2, 3, • • • , k—1 













r -1 








1_ -7y 








< 

1 








1 

7 



— l,o//^Ar — 1,0> 

l^kfi = ^t-1,0 (1- \lk P) 



If the polynomials of degree less than n formed by the first n terms of 

Oq, /3q are input, then the progressive Schur algorithm calculates 
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and 7;. , = 1, • • • , 71 using roughly n - multiplications and additions. 

In contrast with the Szego recursions, Schur’s algorithm provides a class of 
fast Toeplitz solvers based on the Choleski factorization of A/, M = L D L \ with 
L unit lower triangular (see, for example, [19, 22, 3]). In particular, we have the 
following. 

Proposition 2.1. If M = ~ positive definite, then 

Q'q y = — ^y + j, (3qj = /Uy, j = 0,1, . . . , 77—1, are the first n coefficients of 
formal power series a'o,/?o such that = O:q/0q is a Schur function. 

Proposition 2.2. Let M, be as in Proposition 2.1. Then 

1. The Schur parameters obtained by Schur’s algorithm applied to 

are identical with the Schur parameters generated by the Szego re- 
cursions applied to M . 

2. The low'er triangular matrix LD =[’'yx-]yi=o given by 

Tj^k=PkJ-k 

We remark that the connection between Schur functions and Toeplitz matrices is 
provided by the solution of the Carathebdory coefficient problem and the 
correspondence between Carathebdory functions and Schur functions. See. for 
example. [16,1]. Proposition 2.2 follows from by comparison of the resulting 
Szego and Schur recursions. 

Schur’s algorithm can also be formulated in terms of the LFTs T,^ . It is 
shown in [2, 3] that 



T. (r) - 



Cn T Vn ^ 
Vn ”b ^ 



where are polynomials that satisfy the recurrence relations 







Vn—\ 


?n-l 




111 




Co 




0 






In-l 


Vn-l_ 




1 


> 


.^0. 




1 



and |„(X):=X” C„(l/X), (X) : = X" ?/„(l/X). We refer to and >]n- 
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have degree less than n for n >1, as the ?ith Schur polynomials associated with 
the Schur function (})q. The Schur polynomials are generalizations of the Schur 
parameters in the sense that they determine compositions of the elementary linear 
fractional transformations that the Schur parameters determine. 

The Schur polynomials are also important in that they can be used to obtain 
the Szego polynomials. 

Proposition 2.3. Let M, be as in Proposition 2.1, and let be 

the A:th Schur polynomials determined by Then the monic Szego 

polynomials corresponding with M are determined by 

Xki^) ="Vk 



2.3. The Generalized Schur Algorithm. Let 

'^n,k ^n + l°^n+2° ° ^n+k 

be the LFT generated by k steps of Schur’s algorithm applied to 4>„ , and let 
Cn k^Vn k be the Schur polynomials corresponding with k- Note that for any 
> 0 we have <?!>„ = Tj^ k{<i>n-¥k)'i ^bat is, the kth tail of is the (n+^')th tail 
of 6. 

The Generalized Schur Algorithm is based on a doubling procedure for gen- 
erating Ton fi'om . It can be described in general terms as follows. 
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Generalized Schur Algorithm. 

To obtain from 

0. Obtain C„, Vn 

1. Obtain /?(") from , Vn using (f)„ = T~\4 >q), 

2. Obtain Cn,nd/n,n from as were obtained from 

oi'‘\ 

3. Obtain C 2 n^^ 2 n from , C„, using T 2 „ = To „ 



The algorithm is started by performing Step 0 directly by, for example, set- 
ting ^1 = ccqqK^qq and ??] =1. More generally, we can use Schur’s algorithm to 

generate [7y]”° for a small value of Uq, and obtain from the Schur param- 

eters and the recursions (2.1). Now’ suppose Step 0 has been completed for some 
n >?2o- 

^Ve have 






s ^o^ln /^o^n 

■^ = H'Po) = >-■ 

PoVn 



and in [2,3] it is showm that 

^oVn - f^o^n = Tn + l 

PoVn -^bC„ =(5„X" + 0(X" + '). 



We can therefore take 

= (Oo^?n - /^o^n )/ = Wh ~ ^'oC„ )/ X” 

These equations allow us to obtain Q^"^, from Q^*"\ » Vn ui Step 1. 

Step 2 is the doubling step. We obtain „ and ?/„ „ from aj" ) and /?,(") as 
Co.ri = ^ud ?/o,n = ''?n obtained from Q(^"^ and 
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Step 3, the composition of T„ and T„ is given by 

^0,2n Vo,n^n,n ~^^0,nVn,n^ 

^?0,2n ^0,n^n,n d” %,n^ln,n- 

Thus, the Generalized Schur Algorithm consists of various polynomial multi- 
plications and additions performed in a recursive manner. The multiplication of 
polynomials can be efficiently performed using standard FFT techniques. This 
results in an O(nlog“?i) algorithm for computing and where n =2^. The 
Toeplitz system of equations x = b can then be solved by forming Xn from 

by Proposition 2.3, and using the Gohberg-Semencul formula to perform 
the solution phase of the algorithm in 0 (n log?i) arithmetic operations. 

The recursive nature of the algorithm is easily handled in a programming 
language that allows self-referencing subroutines (e.g., PASCAL or APL). In such 
an environment it is easy to obtain a rudimentary implementation of the algo- 
rithm. However, in order to obtain a FORTRAN implementation of the algo- 
rithm we must consider the doubling procedure and resulting storage requirements 
more carefully. 

To contrast the Generalized Schur Algorithm with Schur’s algorithm, the 
coefficients of a^, pj. that are calculated by each algorithm are displayed in Figure 
1 (for no = l). While the classical Schur algorithm generates 7 ^., 

0 < k <n from the number of coefficients of a^., Pf. calculated by the 

Generalized Schur Algorithm is equal to the largest power of two that divides k, 
unless this power of two is less than Uq. Thus, the Generalized Schur Algorithm 
computes only pieces of the Choleski factorization of the Toeplitz matrix. 
Nevertheless, we can obtain every Schur 'parameter; these parameters are often of 
physical and mathematical significance. 



FIGURE 1 



The kth row of squares represents the coefficients of and that are calculated by Schur’s algo- 
rithm and the Generalized Schur Algorithm (A* =0,...,n— 1). The filled squares represent the quanti- 
ties that are calculated during one iteration of each algorithm. In Schur’s algorithm the 
coefficients of Qq, /?o processed one at a time. In the Generalized Schur Algorithm the 

coefficients of Oq are processed in segments, each segment being twice the size of the previous 
segment. 
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3. Numerical Experiments 



We have a working FORTRAN code for the Generalized Schur Algorithm for 
real data based on the procedure described in [4]. Our code implements the Gen- 
eralized Schur Algorithm in real arithmetic, and uses dual split-radix FFT algo- 
rithms to avoid the explicit use of the bit-reversal i)ermutation. The code uses 
6n +0(1) real storage locations, and an additional 3?i/ 2 real storage locations 
for the Toeplitz matrix and the roots of unity (sines and cosines) required by the 
FFTs. The nth monic Szego polynomial Xn squared norm , where 

n =2^, are obtained using (8/3)nlogon + 0(7rlogon) real multiplications and 
(l6/3)n logon + O(7ilogo7i) real additions. We plan to make our code available 
to the scientific community soon. 

In [4] it is shown that can be obtained using less than 8 7i log|n total 

real arithmetic operations. To our knowledge, this is the smallest operation count 
obtained for a superfast Toeplitz solver. Since the Szego recursions require more 
than 2 n~ operations, the Generalized Schur Algorithm requires less computational 
■work for n = 2^ > 256. However, our code does not incorporate the considera- 
tions discussed in [4, Section 5]. The operation count for our code is therefore 
somewhat more than the operation count of [4] in the lower order terms. 

When the Gaussian or Choleski factorization is used to solve the sj'stem of 
equations A x=b, the computed solution x satisfies (.4 +E)x = b with 
l|£ II <£■#•(><) lU It where e is the machine precision and 4>{7i) is a slowly growing 
function of 7i . It follows that the residual r — b — A x — E x and the error 
c — A~^E X satisfy 

IIHI<«'^(-0IUIIII4I. Il4l<c(!<»)lt7|lll^-'llllilt 

Cybenko [9, 10] has shown that 4>[n) is of order 7i- for the Szego recursions (i.e.. 
for the Levinson-Durbin algorithm for the Yule-^\'alker equations). As the 
Choleski decomposition method also yields <^(?! ) of order 7i~. we can conclude that 
the accuracy of the Szego recuisions is comparable with that of Choleski s 
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method. Accurate results can therefore be expected from the Szego recursions 
when the positive definite Toeplitz matrix is well conditioned. 

We now present some numerical experiments to compare the performance of 
the Generalized Schur Algorithm with the Szego recursions. We are pleased to 
report that the results we have obtained are so far favorable. In particular, the 
Generalized Schur Algorithm displays roughly the same stability characteristics 
(e.g., the rates of growth in computed residuals and errors) as the Szego recur- 
sions. The experiments were performed on the VAX 11/750 at Northern Illinois 
University. 

We performed our experiments using the leading principal submatrices of two 
classes of real positive definite Toeplitz matrices + i = o == -^4^+1 • 

We refer to these classes as type A and type B. 

*2 

A {n,0) refers to A/„^j where f-ij = 0^ , — 1<0<1. These matrices have a par- 
ticularly special structure [l4]. The Schur parameters are given by ={—0y- 
For fixed n these matrices become ill-conditioned as 9 tends to ±1. For fixed 9 
the condition numbers tend toward a limiting value as n increases. These are 
numerically banded matrices, where the bandwidth depends on 9 and the smallest 
positive machine number. 

B{n ,9^,92, p) refers to the autocorrelation matrix A/„^j of order /i-fl of a vec- 



(0 <^i < ^2 2)) and rjj is randomly generated in \—p,p\. In other words. 



tor X = [Cy]o^" ^ given by 



Cy = ^y + cos(jw^), 
/fc-i 




where the are n/4 equispaced integral multiples of tt/ 8n in [^ 1 ,^ 2 ] 



Af„^.i = T, where T is the Toeplitz matrix formed by the first n-|-l columns 
of the circulant matrix whose first column is x. These matrices become ill- 
conditioned as p tends to zero, and the ill-conditioning accelerates as ^2 “ 
decreases. 



We first compare the accuracy of the Szego recursions (Sz) and the General- 
ized Schur Algorithm (GS) for computing the monic Szego polynomial x„ and the 
Schur parameters ['7y]". In Tables 1 and 2 we give the square root of the sum of 
the squared errors (the 2-norm of the error) as well as the maximum absolute 
error in the coefficients of Xn and [qfy] ” for the indicated matrices of type A and 
type B. The results of the Szego recursions in double-precision arithmetic are used 
as exact answers for the error calculation. Note that the errors for the two algo- 
rithms seem to grow at roughly the same rate as n increases. In well-conditioned 
problems, the Generalized Schur Algorithm typically results in a somewhat larger 
error than the Szego recursions, but usually this difference is less than an order of 
magnitude. On the other hand, the Generalized Schur Algorithm occasionally 
yielded more accurate results in the ill-conditioned problems of type B. 



TABLE 1 
Phase One Errors 



>1(4096,0.6) 

Szego polynomial errors Schur parameter errors 

n 2-norm max 2- norm max 





GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


4 


0 184e-06 


0 219e-06 


0.119e-06 


0.149e-06 


0 802e-07 


0.112e-06 


0 596e-07 


0 89-1 e- 07 


8 


0 3l5e-06 


0 635e-06 


0.209e-06 


0.328e-06 


0.136e-06 


0 213e-06 


0 782e-07 


0 104e-06 


16 


0 353e-06 


0.738e-06 


0.201e-06 


0.328e-06 


0.181e-06 


0.227e-06 


0.782e-07 


0 104e-06 


32 


0 557e-06 


0 780e-06 


0.209e-06 


0.343e-06 


0.229e-06 


0.227e-06 


0 785e-07 


0 104e-06 


64 


0 719e-06 


0.780e-06 


0.179e-06 


0.343e-06 


0.297e-06 


0.227e-06 


0.988e-07 


0 104 e- 06 


128 


0 983e-06 


0 780e-06 


0.184e-06 


0.343e-06 


0.371e-06 


0 227e-06 


0 988e-07 


0 104e-06 


256 


0.120e-05 


0 780e-06 


0.238e-06 


0.343e-06 


0.436e-06 


0.227e-06 


0 988e-07 


0 104e-06 


512 


0 149e-05 


0.780e-06 


0.358e-06 


0.343e-06 


0.5l2e-06 


0.227e-06 


0 988e-07 


O 

o 

'? 

o 


1024 


0.175e-05 


0.780e-06 


0.477e-06 


0.343e-06 


0.579e-06 


0.227e-06 


0 988e-07 


0 104e-06 


2048 


0.199e-05 


0 780e-06 


0.596e-06 


0.343e-06 


0.639e-06 


0.227e-06 


0.988e-07 


0 104e-06 


4096 


0.225e-05 


0 780e-06 


0.775e-06 


0.343e-06 


0.698e-06 


0.227e-06 


0 988e-07 


0 104e-06 



n 


>1(4096,0 75 ; 

Szego polynomial errors 
2-norm max 


1 

2-1 


Schur parameters errors 
norm max 




GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


4 


0.333e-06 


0.720e-06 


0.238e-06 


0.536e-06 


0.942e-07 


0.321e-06 


0.894e-07 


0 298e-06 


8 


0.300e-05 


0.112e-04 


0.155e-05 


0.572e-05 


0.780e-06 


0.168e-05 


0.745e-06 


0.991e-06 


16 


0.475e-04 


0.276e-04 


0.187e-04 


O.llle-04 


0.473e-05 


0.249e-05 


0.205e-05 


0 991e-06 


32 


0.521e-04 


0.292e-04 


0.197e-04 


0.113e-04 


0.502e-05 


0 252e-05 


0 205e-05 


0 991e-06 


64 


0.59le-04 


0.293e-04 


0.197e-04 


0.113e-04 


0 614e-05 


0 252e-05 


0 205e-05 


0 991e-06 


128 


0 696e-04 


0.293e-04 


0.197e-04 


0 113e-04 


0 736e-05 


0 252e-05 


0.205e-05 


0 991e-06 


256 


0 835e-04 


0.293e-04 


0.197e-04 


0.113e-04 


0 866e-05 


0.252e-05 


0.205e-05 


0 991e-06 


512 


0.986e-04 


0.293e-04 


0.198e-04 


0.113e-04 


0 lOOe-04 


0.252e-05 


0 205e-05 


0 991e-06 


1024 


0.120e-03 


0 293e-04 


0.198e-04 


0.113e-04 


0.118e-04 


0.252e-05 


0 205e-05 


0.991e-06 


2048 


0.140e-03 


0 293e-04 


0.198e-04 


0,113e-04 


0 135e-04 


0.252e-05 


0 205e-05 


0 991e-06 


4096 


0.155e-03 


0 293e-04 


0 199e-04 


0.113e-04 


0.148e-04 


0 252e-05 


0 205e-05 


0 991e-06 



n 


2- 


>1(4096,0.83) 

Szego polynomial errors 

norm max 2-i 


Schur parameters errors 
norm max 




GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


4 


0.449e-05 


0.190e-05 


0.286e-05 


0.143e-05 


0 136e-05 


0 786e-06 


0 131e-05 


0 775e-06 


8 


0.369e-04 


0 173e-03 


0.224e-04 


0.918e-04 


0 420e-05 


0.147e-04 


0 229e-05 


0 115e-04 


16 


0.l80e-02 


0 290e-02 


0 785e-03 


0.114e-02 


0 876e-04 


0.775e-04 


0 643e-04 


0 321e-04 


32 


0.231e-01 


0 636e-02 


0 598e-02 


0.l89e-02 


0.474e-03 


0.112e-03 


0 l57e-03 


0 321e-04 


64 


0.309e-01 


0,655e-02 


0 678e-02 


0.190e-02 


0.638e-03 


0 113e-03 


0 159e-03 


0 32le-04 


128 


0.444e-01 


0 655e-02 


0.678e-02 


0.190e-02 


0 894e-03 


0 113e-03 


0 159e-03 


0 321e-04 


256 


0.547e-01 


0 655e-02 


0.679e-02 


0 190e-02 


0.109e-02 


0 113e-03 


0.159e-03 


0 321e-04 


512 


0.629e-0l 


0 655e-02 


0.680e-02 


0.l90e-02 


0 122e-02 


0 113e-03 


0 lS9e-03 


0 32le-04 


1024 


0.682e-01 


0 655e-02 


0.681e-02 


0 l90e-02 


0 131e-02 


0 1 13e-03 


0 159e-03 


0 321e-04 


2048 


0 725e-01 


0.655e-02 


0 681e-02 


0.190e-02 


0 140e-02 


0 113e-03 


0 159e-03 


0 321e-04 


4096 


0 769e-01 


0 655e-02 


0 682e-02 


0 190C-02 


0 148e-02 


0 113e-03 


0 159e-03 


0 321e-04 



TABLE 2 
Phase One Errors 



j5(4096,20 ',160 ',10.0) 

Szego polynomial errors Schur parameters errors 

n 2-norm max 2- norm max 





GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


4 


0.507e-08 


0.314e-07 


0 301e-08 


0.310e-07 


0.999e-08 


0.999e-08 


0.952e-08 


0.952e-08 


8 


0.880e-07 


0.149e-06 


0.659e-07 


0.838e-07 


0.357e-07 


0.933e-07 


0.242e-07 


0.838e-07 


16 


0.422e-06 


0.850e-06 


0.185e-06 


0.411e-06 


0.134e-06 


0.287e-06 


0.888e-07 


0.158e-06 


32 


0.776e-06 


0.878e-06 


0.294e-06 


0.404e-06 


0.208e-06 


0.340e-06 


0.888e-07 


0.158e-06 


64 


0.186e-05 


0.918e-06 


0.532e-06 


0.387e-06 


0 499e-06 


0.413e-06 


0.151e-06 


0,158e-06 


128 


0.256e-05 


O.llle-05 


0.841e-06 


0.504e-06 


0.587e-06 


0.434e-06 


0.151e-06 


0.158e-06 


256 


0.399e-05 


0 138e-05 


0.138e-05 


0.692e-06 


0.683e-06 


0.477e-06 


0.151e-06 


0 158e-06 


512 


0.670e-05 


0.181e-05 


0.238e-05 


0.734e-06 


0.781e-06 


0.540e-06 


0.151e-06 


0.158e-06 


1024 


0.115e-04 


0 246e-05 


0.441e-05 


0.955e-06 


0.104e-05 


0.656e-06 


0.151e-06 


0.158e-06 


2048 


0.205e-04 


0.344e-05 


0.811e-05 


0.966e-06 


0.139e-05 


0.919e-06 


0.151e-06 


0.158e-06 


4096 


0.854e-04 


0.270e-04 


0.297e-04 


0.119e-04 


0.547e-04 


0.168e-04 


0.250e-04 


0.616e-05 







B(4096,20 *,160 
Szego polynomial errors 


*,0.5) 


Schur parameters errors 




n 


2- norm 


max 


2- norm 


max 




GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


4 


0.346e-07 


0.587e-07 


0.257e-07 


0.515e-07 


0.328e-07 


0.570e-07 


0.239e-07 


0.515e-07 


8 


0.945e-07 


0.474e-06 


0.596e-07 


0.232e-06 


0.625e-07 


0.261e-06 


0.527e-07 


0.232e-06 


16 


0.807e-05 


0.377e-04 


0 284e-05 


0.144e-04 


0.116e-05 


0.527e-05 


0.873e-06 


0.430e-05 


32 


O.llle-02 


0.480e-02 


0.406e-03 


0.195e-02 


0.467e-04 


0.205e-03 


0 204e-04 


0.937 e- 04 


64 


0.436e-02 


0.499e-02 


0.172e-02 


0.203e-02 


0.148e-03 


0.312e-03 


0.588e-04 


0.999e-04 


128 


0.501e-02 


0.524e-02 


0.156e-02 


0.187e-02 


0.185e-03 


0.362e-03 


0.588e-04 


0.999e-04 


256 


0.666e-02 


0.848e-02 


0.153e-02 


0.163e-02 


0 238e-03 


0.490e-03 


0. 588c-04 


0.999e-04 


512 


0.796e-02 


0.182e-01 


0.158e-02 


0.170e-02 


0.280e-03 


0.809e-03 


0.588e-04 


0.999e-04 


1024 


0.139e-01 


0.316e-01 


0 189e-02 


0.231e-02 


0.428e-03 


0.140e-02 


0.588e-04 


0.999e-04 


2048 


0.461e-01 


0.659e-01 


0.303e-02 


0.353e-02 


0.145e-02 


0.311e-02 


0.116€-03 


0.173e-03 


4096 


0.506e-01 


0.413e+00 


0.131e-01 


0.104e+00 


0.376e-01 


0.221e-h00 


0.l32e-01 


0 830e-01 






5(4096,50 ‘,100 
Szego polynomial errors 


',1.0) 


Schur parameters errors 




n 


2-norm 


max 


2- norm 


max 




GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


4 


0.393e-05 


0.409e-05 


0.238e-05 


0.262e-05 


0.169e-05 


0.177e-05 


0.161e-05 


0.173e-05 


8 


0.490e-04 


0.303e-04 


0.312e-04 


0.277e-04 


0.270e-04 


0.255e-04 


0.206e-04 


0.187e-04 


16 


0, 124e-03 


0.307e-03 


0.830e-04 


0.154e-03 


0.585e-04 


0.120e-03 


0.427e-04 


0.631e-04 


32 


0.319e-03 


0.500e-03 


0.140e-03 


0.170e-03 


0.129e-03 


0.224e-03 


0.480e-04 


0.757e-04 


64 


0.411e-03 


0.647e-03 


0.120e-03 


0.165e-03 


0.177e-03 


0.333e-03 


0.480e-04 


0.757e-04 


128 


0.644e-03 


0.890e-03 


0.141e-03 


0.197e-03 


0.245e-03 


0.479e-03 


0.634e-04 


0.757e-04 


256 


0.952e-03 


0.129e-02 


0.165e-03 


0.246e-03 


0.346e-03 


0.676e-03 


0.634e-04 


0.783e-04 


512 


0.152e-02 


0.231e-02 


0.196e-03 


0.309e-03 


0.530e-03 


0.104e-02 


0.634e-04 


0.989e-04 


1024 


0.260e-02 


0 325e-02 


0.231e-03 


0.357e-03 


0.904e-03 


0.155e-02 


0.726e-04 


0.114e-03 


2048 


0.303e-02 


0.439e-02 


0.254e-03 


0.396e-03 


0.122e-02 


0.217e-02 


0,726e-04 


0.114e-03 


4096 


0.395e-02 


0.641e-02 


0.282e-03 


0.436e-03 


0.173e-02 


0.312e-02 


0.769e-04 


0.114e-03 
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Oui’ second accuracy experiment compares the Szego recursions and the Gen- 
eralized Schur Algorithm for the solution of a Toeplitz system of equations. The 
Gohberg-Semencul formula is used for the solution phase for both algorithms. 

Let ||M||i and ||Af||2 denote the norm of the matrix M that is induced by 
the vector 1-norm and vector 2-norm, respectively. Then 

max ||M e,- II2 < ||M||o < ||A/||i. (3.1) 

1 < y < n 

Let A4 = and define R : = R D~'^ so that M~^ = R R* . It is shown in [2] 
that 

IlfieJII <l|V->|l2<n ||A ej||. (3.2) 

This gives bounds on the 2-norm of M~^ in terms of the last column of R (i.e., in 
terms of the coefficients of the normalized Szego polynomial )• 

After computing a solution x of the positive definite Toeplitz system, the 
residual r = b — M x \s calculated in double precision, and the system AJ e = r 
is solved (in 0(?! Ig n) operations) to obtain an estimate e of the error. Finally. 

Ilrll 

the upper and lower bounds for the two ratios = -7^ — — ajuj 

llAUIIIill 

■ 4 ^, = — — , — are obtained using (3.1) and (3.2). In order for an algo- 

IlMlI \\M-^\\ ll-^ll 

rithm to be stable, these ratios should increase slowly with n. In Tables 3 and -4 
we display the computed lower bounds on computed relative residu- 

als and relative errors, and upper bounds on the ratios and ipf These results 
were obtained using b = [y-fl]”^o- Similar results were obtained when the com- 
ponents of b were set equal to one and when they were randomly generated. The 
Szego recursions were stopped at ?i = 2 . 

Note that the sizes of the residuals and errors for each algorithm are roughly 
the same, as are the rates of gi'owth in these quantities. The upper bounds on i/v 
and ipg are slowly growing for both algorithms. The errors and residuals for the 



n 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

2048 

4096 

8192 

16384 

32768 

n 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

2048 

4096 

8192 

16384 

32768 

n 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

2048 

4096 

8192 

16384 

32768 



TABLE 3 

Toeplitz Solver Residuals and Errors 



II2 


Hr II 

imr 




32 .— — 






V'', 






GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


0.449e+01 


0.3T6e-06 


0.580e-06 


0.353e-06 


0.483e-06 


0.376e-06 


0.580^-06 


0.594e-07 


0.812e-07 


0.493e+01 


0.22le-06 


0.6l5e-06 


0.253e-06 


0.631e-06 


0.266e-06 


0.740e-06 


0.388e-07 


0.967e-07 


0.495e+01 


0.3S3e-08 


0.444e-06 


0.594e-06 


0.898e-06 


0.545e-06 


0.63le-06 


0.906e-07 


0.137^06 


0.495e+0l 


0.2SOe-06 


0.457e-06 


0.540e-06 


0.677e-06 


0.449e-06 


0.731e-06 


0.825e-07 


0.l03e-06 


0.495e+01 


0.725e-06 


0.117e-06 


0.910e-06 


0.711e-06 


0.125e-05 


0.202e-06 


0.139e-06 


0.109^06 


0.495e+01 


0.567e-06 


0.296e-06 


0.l20e-05 


0.107e-05 


0.102e-05 


0.529e-06 


0.184e-06 


0.163e-06 


0.495e+01 


0.l05e-05 


0.205e-06 


0.l7le-05 


0.140e-05 


0.192e-05 


0.375e-06 


0.26le-06 


0.214e-06 


O.495e+0l 


0.l50e-05 


0.327e-06 


0.161e-05 


0.697e-06 


0.278e-05 


0.604e-06 


0.246e-06 


0.106e^06 


0.495e4-0l 


0.206e-05 


0.399e-06 


0.222e-05 


0.893e-06 


0.383e-05 


0.743e-06 


0.338e-06 


0.136-06 


0.495e+01 


0.2S8e-05 


0.446e-06 


0.299e-05 


0.934e-06 


0.538e-05 


0.833e-06 


0.457e-06 


0.143^06 


0.495e-f01 


0.410e-05 


0.501e-06 


0.4l8e-05 


0.944e-06 


0.767e-05 


0.937 e-06 


0.638e-06 


0.144e-06 


0.495e+01 


0.595e-05 


0.683e-06 


0.600e-05 


0.109e-05 


O.llle-04 


0.128e-05 


0.916e-06 


0.167-06 


0.495e+01 


0.571e-05 




0.590e-05 




0.107e-04 




0.90le-06 




0.495e4-0l 


0.723e-05 




0.738e-05 




0.135e-04 




0.113^05 





>l(2'^0.75) 



uai 1I2 


Hr II 
W 


II 

II 


HI 
i II 






V'-, 


r 




GS 


Sz 


GS 


Sz 


GS 


Sz 


GS 


Sz 


0.279e+02 


0.703e-06 


0.122e-05 


0.434e-06 


0.l69e-05 


0.321e-06 


0.557e-06 


0.102e-07 


0.397e-07 


0.647e+02 


0.704e-06 


0.591e-06 


0.278-05 


0.514e-05 


0.316e-06 


0.265e-06 


0.281e-07 


0.520e-07 


0.799e+02 


0.170e-05 


0.163e-05 


0.905e-05 


0.115e-04 


0.967e-06 


0.923e-06 


0.741e-07 


0.942e-07 


0.803e+02 


0.123e-05 


0.425e-05 


0.140e-04 


0.117e-04 


0.959e-06 


0.33le-05 


0.114e-06 


0.954e-07 


0.803e+02 


0.386e-05 


0.285e-05 


0.161e-04 


0.142e-04 


0.402e-05 


0.296e-05 


0.131e-06 


0.116e-05 


0.803e+02 


0.5l3e-05 


0.293e-05 


0.289e-04 


0.260e-04 


0.681e-05 


0.389e-05 


0.235^^06 


0.212e06 


0.803e+02 


0.922e-05 


0.255e-05 


0.520e-04 


0.504e-04 


0.148e-04 


0.408e-05 


0.424^06 


0.411e-06 


0.803e+02 


0.384e-05 


0.286e-05 


0.267e-04 


0.207e-04 


0.699e-05 


0.52le-05 


0.2l7e-06 


0.169e-05 


0.803e+02 


0.6l0e-05 


0.300e-05 


0.372-04 


0.315e-04 


0.120e-04 


0.590e-05 


0.303e-06 


0.257e-06 


0.803e+02 


0.89le-05 


0.288e-05 


0.363-04 


0.355e-04 


0.183e-04 


0.593e-05 


0.296e-06 


0.289e-06 


0.803e+02 


0.951e-05 


0.322e-05 


0.317e-04 


0.303e-04 


0.200e*04 


0.679e*05 


0.258e-06 


0.246-06 


0.803e+02 


0.692e-05 


0.334e-05 


0.338e-04 


0.341e-04 


0.148e-04 


0.714e-05 


0.275e-06 


0.278e<06 


0.803e+02 


0.l07e-04 




0.744^04 




0.231e-04 




0.606- 06 




0.803e+02 


0.1o9e-04 




0.738e-04 




0.344e-04 




0.601e-06 
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0.72le+02 


0.161e-05 


0.265e-05 


0.741e-06 


0. 129-05 


0.448-06 


0.739-06 


0. 636-08 


0.110-07 


0.320e+03 


0.461e-05 


0.546e-05 


0.661e-05 


0.213e-05 


0.932-06 


0.110-05 


0.127-07 


0. 408-08 


0.610e+03 


0.462e-05 


0.l24e-04 


0.327e-04 


0.826-05 


0.985-06 


0.265-05 


0.329-07 


0.831-08 


0.642e+03 


0.764e-05 


0.213e-04 


0.137e-03 


0.190-04 


0.225-05 


0.628-05 


0.131-06 


0.182-07 


0.642e+03 


0.1l2e-04 


0.192e-04 


0.162-03 


0.400-04 


0.468-05 


0.800-05 


0 155-06 


0.383- 07 


0.642e+03 


0.427e-05 


0.154e*04 


0. 190-03 


0.123-03 


0.249-05 


0.900-05 


0.181-06 


0.118-06 


0.642e+03 


0.l26e-04 


0.172e-04 


0.279-03 


0.291-03 


0,101-04 


0.138-04 


0.267-06 


0.278-06 


0.642e+03 


0.113e-04 


0.174e-04 


0.264-03 


0.129-03 


0.121-04 


0.186-04 


0.253-06 


0.123-06 


0.642e+03 


0.876e-05 


0.172e-04 


0.363-03 


0.264-03 


0.120-04 


0.237-04 


0.347-06 


0.253-06 


0.642e+03 


0.486e-05 


0.175e-04 


0.392-03 


0.301-03 


0.812-05 


0.292-04 


0.375-06 


0.288- 06 


0.642e+03 


0.817e-05 


0.150e-04 


0.367-03 


0.295-03 


0.156-04 


0.288-04 


0.351-06 


0.282-05 


0.642e+03 


0.1l6e-04 


0.178e-04 


0.421-03 


0.372-03 


0.241-04 


0.371-04 


0.403-06 


0.356-06 


0.642e+03 


0.221e-04 




0.891-03 




0.482-04 




0.852-06 




0.642e+03 


0.204-04 




0.906-03 




0.458-04 




0.867-06 





n 

4 

8 

16 

32 

54 

128 

256 

512 

1024 

2048 

4096 

8192 

16384 

32768 

n 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

2048 

4096 

8192 

16384 

32768 

n 

4 

8 

16 

32 

64 

128 

256 

512 

1024 

2048 

4096 

8192 

16384 

32768 



TABLE 4 

Toeplitz Solver Residuals and Errors 
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0.133e+01 


0.824e-07 


0.824e-07 


0.760^07 


0.760e-07 


0.545^07 


0.545^07 


0.538e-07 


0.538e-07 


0.178e+01 


0.211e-06 


0.235^06 


0.l36e-06 


0.249e-06 


0.828e-07 


0.922e-07 


0.701e-07 


0.129e-06 


0.l92e+01 


0.372^06 


0.221e06 


0.112e-06 


0.170e06 


0.763e-07 


0.453^07 


0.529e-07 


0.806e-07 


0.200e+01 


0.5l5e-06 


0.305e06 


0.468e-06 


0.173e-06 


0.889e-07 


0.526e07 


0.212e-06 


0.787e-07 


0.203e+01 


0.408^06 


0.248e-06 


0.198e-06 


0.878e-07 


0.657^07 


0.399e-07 


0.883e-07 


0.391 e-07 


0.205e+01 


0.538e-06 


0.389e-06 


0.270e-06 


0.224e-06 


0.840^07 


0.607 e-07 


0.119e-06 


0.988^07 


0.206e+01 


0.101^05 


0.736^06 


0.904e-06 


0.650e-06 


0.156^06 


0.114^06 


0.396e-06 


0.285e-06 


0.206e+01 


0.395e-05 


0.771e-06 


0.393e-05 


0.677e-06 


0.612e-06 


0.120e-06 


0.172e-05 


0.296^06 


0.207e+01 


0.115e-04 


0.1l7e-05 


0.1l5e-04 


0.110^05 


0.175e-05 


0.177e-06 


0.501e-05 


0.482e-06 


0.208e+01 


0.228e-04 


0.514e-05 


0,229e-04 


0.514e-05 


0.339e-05 


0.763^06 


0.992e-05 


0.223e-05 


0.209e+01 


0.50l€-04 


0.164e-04 


0.501e-04 


0.164e-04 


0.723e-05 


0.237 e-05 


0.216«^04 


0.708e-05 


0.2l2e+01 


0.102^03 


0.375e-04 


0.102e-03 


0.375e-04 


0.144e-04 


0.532^05 


0.432e-04 


0.l60e-04 


0.2l7e+01 


0.196e-03 




0.197 e- 03 


O.OOOe+00 


0.265e-04 


O.OOOe+00 


0.814e-04 


O.OOOe+00 


0.614e+01 


0.469e-03 




0.470e-03 
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0.399e-04 
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0.520e-04 


O.OOOe+00 
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0.158e+01 


0.847e-07 


0.942e-07 


0.964e-07 


0.119e-06 


0.484e-07 


0.539e-07 


0.560e-07 


0.689^07 


0.382e+01 


0.304^06 


0.262^06 


0.192^06 


0.171e-06 


0.732e-07 


0.631 e- 07 


0.448e-07 


0.398e-07 


0.440e+02 


0.486e-05 


0.314^05 


0.338e-05 


0.559e-06 


0.137e-06 


0.886e-07 


0.676e-07 


0.1l2e-07 


0.974e+02 


0.273e-04 


0.272e-04 


0.201^04 


0.112e-04 


0.696e-07 


0.693e-07 


0.l80e-06 


O.lOOe-05 


0.112e+03 


0.391e-04 


0.285e-04 
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0.8l0e-07 


0.590e-07 


0.207e-06 


0.125e-06 


0.118e+03 


0.423e-04 


0.363e-U4 


0.299^04 


0.194e-04 


0.800e-07 


0.687e-07 


0.220e-06 


0.143e-06 


0.121e+03 


0.600e-04 


0.426e-04 


0.471e-04 


0.20le-04 


0.l09e-06 


0.777e-07 


0.337e-06 


0.144e-06 


0.123e+03 


0.729e-04 


0.457e-04 


0.662e-04 


0.196e-04 


0.132e-06 


0.827e07 


0.468e-06 


0.138e-06 


0.l24e+03 


0.914e-04 


0.484e-04 


0.830e-04 


0.189e-04 


0.161e-06 


0.855e-07 


0.582e-06 


0.133e-06 


0.124e-h03 


0.933^04 


0.456e-04 


0.833e-04 


0.144e-04 


0.160e-06 


0.783e-07 


0.580e-06 


O.lOle-OG 


0.126e+03 


0.675^04 


0.472e-04 


0.500e-04 


0.705e-05 


0.113e-06 


0.787e-07 


0.345e-06 


0.487e-07 


0.128e+03 


0.607^04 


0.753e-04 


0.235e-04 


0.530e-04 
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0.276e+03 


0.223e-04 


0.188e-04 


0.513e-05 


0.109e-04 


0.717e-07 


0.606^07 


0.113^07 


0.239e-07 


0.771e+04 


0.654e-03 


0.147e-02 


0.273e-03 


0.125e-03 


0.482e-07 


0.109e-06 


0.183e-07 


0.840^08 


0.819e+04 


0.720e-03 


0.148e-02 


0.332e-03 


0.238e-03 


0.342e-07 


0.704e-07 


0.200^07 


0.143e-07 


0.872e+04 


0.118e-02 


0.161e-02 


0.757e-03 


0.662e-03 


0.505e-07 


0.687e-07 


0.417e-07 


0.364e-07 


0.926e+04 


0.222e-02 


0.179e-02 


0.141e-02 


0.986e-03 


0.905e-07 


0.728e-07 


0.722e-07 


0.504 e-07 


0.939e+04 


0.33le-02 


0.208e-02 


0.270e-02 


0.124e-02 


0.132e-06 


0.83le-07 


0.136e-06 


0.619e-07 


0.953e+04 


0.419e-02 


0.209e^02 


0.394e-02 


0.122e-02 


0.165e-06 


0.821e-07 


0.194e-06 


0.600e-07 


0.959e+04 


0.584e-02 


0.196e-02 


0.554e-02 


0.994e-03 


0.226e-06 


0.755e-07 


0.271e-06 


0.483^07 


0.963e+04 


0.669e-02 


0.201^02 


0.645(^02 


0.888e-03 


0.254e-06 


0.759e-07 


0.313^06 


0.430e-07 


0.969e+04 


0.749e-02 


0.202e-02 


0.720e-02 


0.826e-03 


0.283e-06 


0.759e-07 


0.348e-06 


0.398e-07 


0.978e+04 


0.825e-02 


0.l94e-02 


0.785e-02 


0.739e-03 


0.318e-06 


0.743e-07 


0.375^06 


0.352^07 


O.lOle+05 


0.83le-02 


0.l78e-02 


0.117e-01 


0.715e-03 


0.323e-06 


0.684^07 


0.544e-06 


0.336e07 


0.103e+05 


0.802e-02 




0.789e-02 




0.333e-06 




0.358^06 




0.116e+05 


0.804e-02 




0.277e-01 




0.400e-06 




0.nie05 
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Generalized Schur Algorithm were typically larger than those of Szego recursions 
in the well-conditioned problems, while the Generalized Schur Algorithm some- 
times performed slightly better in the relatively ill-conditioned problems. These 
differences were rarely more than an order of magnitude. 

Thus, in the examples we investigated the accuracy of the Generalized Schur 
Algorithm closely followed that of Szego recursions. Of course, while we did not 
encounter any problems where the accuracy of the two algorithms was greatly 
different, more experimentation is needed to obtain a better confidence in the 
results of the Generalized Schur Algorithm. 
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Table 5 shows average CPU times required by the Generalized Schur Algo- 
rithm and Szego recursions for the calculation of the nth Szego polynomial and 
the n Schur parametei-s. Note that the Szego recursions would have used about 36 
CPU houre for n =2^®, compared with about 38 minutes for the Generalized 
Schur Algorithm. The dramatic relative efficiency of the Generalized Schur Algo- 
rithm for large problems and the potential reliability of the algorithm as indicated 
in the above examples lead us to believe the algorithm will become a popular 
method for solving huge positive definite Toeplitz systems of equations. 





TABLE 5 




Timing 


comparison (CPU seconds) 


n 


GS 


Sz 


256 


3.8 


2.4 


512 


8.1 


8.6 


1024 


17.8 


32.2 


2048 


40.0 


126.0 


4096 


90.7 


505.8 


8192 


200.8 


2030.9 


16384 


441.2 


8068.3 


32768 


1001.6 


32222.3 


65536 


2257.4 


— 
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